

%% load estimation results
load('multiple_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z specif l0 nmissing


%% load everything of interest
load('../../baseline_model/sim_eq_round2.mat')
N=length(nmissing);
intervEq=intervEq(:,nmissing);
U0=U0(:,:,nmissing);
V_R=V_R(:,:,nmissing);
V_US=V_US(:,:,nmissing);
V_UK=V_UK(:,:,nmissing);
V_France=V_France(:,:,nmissing);
V_Russia=V_Russia(:,:,nmissing);
V_China=V_China(:,:,nmissing);
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% standard errors via Hessian
InvSigmaH = HessPhat(thetaStar,1,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
seH = sqrt(diag(inv(InvSigmaH)));
resultsH = printResults([thetaStar,seH,2*normcdf(abs(thetaStar./seH),'upper')],reb,mpow,specif);
disp(' ')
disp('for Tables E5-E6:')
disp(resultsH)


%%
clearvars -except thetaStar InvSigmaH seH resultsH
save('inference_multiple.mat')



